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ABSTRACT 

In this work, the gas infall and the formation of outflows around low and high mass proto- 
stars are investigated. A radial self-similar approach to model the transit of the molecular gas 
around the central object is employed. We include gravitational and radiative fields to produce 
heated pressure-driven outflows with magneto-centrifugal acceleration and collimation. Outflow 
solutions with negligible or vanishing magnetic field are reported. They indicate that thermody- 
namics is a sufficient engine to generate an outflow. The magnetized solutions show dynamically 
significant differences in the axial region, precisely where the radial velocity and collimation are 
the largest. They compare quantitatively well with observations. The influence of the opacity on 
the transit solutions is also studied. It is found that, when dust is not the dominant coolant, such 
as in the primordial universe, mass infall rates have substantial larger values in the equatorial 
region. This suggests that star forming in a dust-free environment should be able to accrete 
much more mass and become more massive than present day protostars. It is also suggested 
that molecular outflows may be dominated by the global transit of material around the protostar 
during the very early stages of star formation, especially in the case of massive or dust-free star 
formation. 

Subject headings: ISM: jets and outflows - stars: formation - methods: analytical 



1. Introduction 

The understanding of star formation has changed 
dramatically during the last two decades thanks 
to the development of both instrumental and nu- 
merical techniques. However, while the picture 
of the main stages of low mass star formation is 
becoming clearer, the situation for massive ob- 
jects remains problematic. In the present paper, 
we address the problem of the origin of outflows 
and infall around protostars over a wide range of 
masses. 

At present, the stages of low mass protostar 
evolution are empirically divided in four classes 
(see for example Andre et al. (2000) for a review) 
corresponding to the evolutionary sequence from 
the initial collapse to a new star. "Class 0" proto- 
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stars correspond to the earliest stage of star forma- 
tion. Observationally, they appear deeply embed- 
ded in a circumstellar dusty envelope (detected in 
sub-millimeter wavelengths) that is more massive 
than the central stellar mass. Furthermore, they 
show powerful bipolar ejections of material in the 
form of collimated CO outflows which distinguish 
them from the pre-stellar phase of star formation 
(a gravitationally bound core within a molecular 
cloud) . 

The next stage of the star formation process 
corresponds to "Class I" objects which are typi- 
cally I0 5 year old and are characterized by a pos- 
itive slope of their spectral energy distributions 
(SEDs) between 2.2 and 10-25 /im. They are still 
surrounded by a diffuse circumstellar envelope but 
are also surrounded by a disk from which the ac- 
cretion onto the central object takes place. At 
this stage, the mass contained in the envelope is 
smaller than the mass of the central protostar. 
Like Class objects, they also produce bipolar 
outflows (optical jets and molecular outflows) but 
these are less (~ I order of magnitude) powerful 



1 



than those observed for Class protostars (Bon- 
temps et al. 1996). 

Class II and III objects (resp. classical and 
weak-line T Tauri stars) are the final two stages in 
evolutionary sequence of protostars. These corre- 
spond to pre-main sequence stars surrounded by 
an accretion disk (optically thick (Classical) or 
thin depending on the degree of evolution) but 
which do not have a circumstellar envelope, as op- 
posed to Class or 1 objects. 

Formation of massive stars (M > 10M Q ), is 
less well understood but two main theories have 
emerged: i) the coalescence scenario in which mas- 
sive stars are formed by the merging of interme- 
diate mass objects (Bonnell et al. 1998) and ii) 
the accretion scenario, originally applied to mas- 
sive star by Beech & Mitalas (1994), where a 
massive star forms thanks to large accretion rates 
(1CT 4 -1CT 2 M©) onto the central object (Norbcrg 
& Maeder 2000; Behrend & Maeder 2001). The 
main difficulty of the merging scenario arises when 
one considers the high stellar densities needed for 
it to be efficient. On the other hand, it is difficult 
to accrete gas onto a very luminous star so that 
the upper mass of the accretion scenario is rapidly 
reached. 

Whatever the mass of the central object is, it 
appears that the accretion of material is accompa- 
nied by strong bipolar outflows. During the last 
two decades, starting with the first outflow obser- 
vation from a forming star by Sncll et al. (1980), 
the number of outflow observations has increased 
to the point that this phenomenon is now widely 
believed to affect every forming star. Bipolar out- 
flows allow the forming star to transport any ex- 
cess of angular momentum. They can be classified 
as either atomic jets or molecular outflows, accord- 
ing to their properties. 

Jets are fast (~ 100 - 300 km s" 1 ), well- 
collimated (opening angle < 10°), and mostly con- 
stituted of atomic material. Many models and 
simulations attempt to describe how they can be 
launched from the magnetized accretion disk of the 
protostar (Ferreira 1997; Ouyed & Pudritz 1997), 
and their magneto-centrifugal origin is commonly 
accepted. However, the nature of the exact mech- 
anism, be it a disk wind (Blandford & Payne 1982) 
or the asymptotic collimation of an A-wind (Shu 
ct al. 1995, 2000), is still debated. 



Molecular outflows (often traced by CO and 
H 2 molecules) are more massive, slower (~ tens 
of km s -1 ) and less collimated than jets. The 
driving mechanism of these outflows is still open 
to discussion but the jet-driven bow shock model 
(Raga & Cabrit 1993; Chcrnin ct al. 1994; Downes 
& Ray 1999; Ostrikcr ct al. 2001) and the wind- 
driven model (Shu 1991; Shu ct al. 2000) arc the 
two main theories that have emerged. In the jet- 
driven model, the bow shock surface created at the 
head of the jet transfers momentum to the ambi- 
ent medium thus producing a thin shell that is 
identified as the molecular outflow. In the second 
approach, a wide-angle wind creates two bipolar 
wind-blown bubbles that sweep up the ambient 
material (then referred to as the molecular out- 
flow). 

For low and intermediate mass star formation, 
the morphologies and kinematics of the observed 
outflows are generally explained by one of these 
two models: some outflows present jet-driven fea- 
tures (e.g. HH212, Orion S outflow, Rodriguez- 
Franco et al. (1999)) whereas others (e.g. VLA 
0548) show wind-driven signatures (Lcc ct al. 
2000, 2001). However, when massive star for- 
mation is considered, the observations cannot be 
matched by any of these two models (Churchwcll 
1997, 2000). The formation of massive stars gives 
rise to molecular outflows showing different prop- 
erties than the ones recorded in low and interme- 
diate mass star formation: they are very massive 
(~ tens of solar masses and, in some cases, more 
massive than the central stars that presumably 
drive them), poorly collimated and also faster than 
in the case of low mass star formation. Accord- 
ing to Churchwell (1997, 2000), jets cannot entrain 
more than a few solar masses and certainly not the 
tens of solar masses observed in massive outflows. 
Furthermore, the wide opening angles observed 
are hardly explained by the entrainment from a 
collimated jet. On the other hand, if entrainment 
by a wide-angle wind allows a large opening an- 
gle for the outflow, it cannot match the velocities 
reached by massive outflows. 

In this paper, we present a general radial self- 
similar model for the flows surrounding young 
forming stars, applicable in both low and high 
mass star formation. The paper is laid out as fol- 
lows: in Sect. 2, we describe the model and its 
equations and show corresponding hydrodynami- 
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cal solutions in Sect. 3; some magnetized solutions 
are then presented, together with a study of the 
influence of changes in the opacity, in Sect. 4; we 
give the main properties of the model in Sect. 5 
and discuss them in Sect. 6 before concluding. 

2. Theoretical basis of Transit Models 

The present models are based on ideal MHD 
with the assumption of radial self-similarity and a 
simplified treatment of radiation (Fiege & Henrik- 
scn f996; Lery 2003). 

Spherical coordinates are used: the system is 
centered on the protostar, the x — y plane corre- 
sponds to the equatorial plane and the poloidal 
angle 9 is taken relative to the rotational axis of 
the central object. 

2.1. Fluid equations 

We describe the gas with the usual macroscopic 
quantities: density p, velocity v, pressure p and 
temperature T. The size of the system is such 
that the conditions for ideal MHD are fulfilled. 
We further assume steady-state (d/dt = 0), axial 
symmetry around the rotational axis (d/d(j) = 0). 
The standard ideal MHD equations then reduce 
to: 

V • (pv) = (1) 



-Vp + V$+(v- V)v= VxBxB (2) 

p Airp 

VB = 0. (3) 

The magnetic Reynolds number being much 
greater than unity, the diffusion term in the induc- 
tion equation can be neglected, and the magnetic 
field lines are frozen in the plasma. The induction 
equation becomes 

V x (v x B) = (4) 

which directly implies 

v x B = V* (5) 

where >J> is the electric potential (E = — Vf). A 
non-zero toroidal electric field could convect the 
magnetic poloidal field lines into a sink on the axis 



of the system. To ensure a zero toroidal electric 
field, the only requirement from Eq. (5) is v p oc B p 
(Chan & Hcnrikscn 1980; Henriksen f996), where 
the subscript p denotes the poloidal component of 
the field (v = v p + and B = B p + B^). 

We also consider the equation of state for ideal 
gas, i.e. p = nkT. 

Note also that the purely hydrodynamical prob- 
lem will be studied by using Eq.(f ) and setting the 
RHS in Eq.(2) to zero. 

2.2. Self-similar treatment 

2.2.1. Self- similarity - general statements 

A phenomenon is called self-similar if the spa- 
tial (or temporal) distributions of its properties at 
various different times (or locations) can be ob- 
tained from one another by a similarity transfor- 
mation. The investigation of the full phenomenon 
can then be reduced to the study of the properties 
of the system for only a specific time (or loca- 
tion). If the origin of time can be chosen arbitrar- 
ily, the scales of length and mass are also arbitrary, 
and the system is 'scale-free'. The system of par- 
tial differential equations (PDE) that describe the 
problem is transformed to a set of ordinary differ- 
ential equations (ODE), which drastically simpli- 
fies the investigation. 

Self-similar models turn out not only to de- 
scribe the behavior of physical systems under 
some special conditions, but also describe the 
intermediate-asymptotic behavior of solutions to 
wider classes of problems in the range where these 
solutions no longer depend on the details of the 
initial and/or boundary conditions, yet the sys- 
tem is still far from being in an ultimate equilib- 
rium state (Barenblatt & Zel'dovich 1972). For 
an extensive review of self-similarity, the reader is 
referred to Barenblatt (1996). In present case of 
star formation, we express each variable using the 
following self-similar form 

L(r,6)=L x (j-^J xl(6). (6) 

The constant Lq has the dimensions of the variable 
and is directly dependent on the parameters of the 
problem such as the mass or the luminosity of the 
central object. The free parameters are the self- 
similar indices a and the fiducial scale r . 
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2.2.2. Self- similar fluid quantities 

Several observations of the protostcllar environ- 
ment favor our self-similar approach since they 
show that radial density profiles in regions of iso- 
lated star formation are well represented by a 
power-law p oc r~ p with p ~ 0.5 — 2 (e.g. L1527 
has been found to have p in the range of 1.5-2 
Ladd ct al. (1991)). Indeed, in some cases, den- 
sity profiles are reproduced by most hydrostati- 
cally supported (p = 2) or free infalling (p = 1.5) 
cloud core, e.g. Shu (1977). In some other cases, 
the slope has been found to be much shallower 
than the ones predicted by the standard models. 
Density profiles as shallow as 0.5-0.9 can be in- 
terpreted as the presence of a magnetic (Barsony 
& Chandler 1993) or rotational support (Chandler 
et al. 1998) of the core. 

We now apply our assumptions to the ideal 
MHD equations. Dimensional analysis allows us 
to determine the index of each variable as a func- 
tion of a single free self-similar parameter, a. The 
systems then reads 



p(r,0) 



p(r,6) = 



M 



GM Z 



ro 



P{6) 



T(r,6) 



/j,rriH GM 



< m (GM\ 



Q(6) 



- u(0) 



(7) 

(8) 
(9) 

(10) 



B r ,e,<p( r , 6) 



GM< 



ro 



u r fl,4>{0) 
Dp,p,4>{0) 



(11) 

In these expressions, G is the gravitational con- 
stant, ks the Boltzmann constant, m H the mass 
of the proton, p, the mean molecular weight (we 
use p, = 2), and M the mass of the central ob- 
ject. The fiducial scale r is a radius of refer- 
ence. The self-similar assumption is valid only 
above this radius (see Sect. 2.3.2). Note also that 
in Eq. (11), the subscript p stands for "poloidaF. 
The self-similar index a is a free parameter that 
lies between —1/2 < a < 1/4. In particular, 
a = —1/2 yields to pure radial accretion. This 
issue is discussed in Ficgc & Hcnrikscn (1996) 



and the behavior of the solutions with a can be 
found in Lery et al. (1999). Note that the hydro- 
static (p = —2a + 1/2 = 2) and the free infalling 
(p = 1.5) cases cannot be treated by the present 
model. 

For the remainder of this paper, we choose 
a = —0.2 in most cases. From Eq.(7), we infer 
p oc r~ - 9 , which corresponds to a shallow den- 
sity profiles where rotation and magnetic field can 
drive substantial effects. 

In a nearly study on self-similar transit mod- 
els (Fiege & Henriksen 1996), the same propor- 
tionality relationship between each component of 
the velocity and magnetic field was taken. All the 
components of the electric field Eq.(5) were then 
equal to zero. Such a configuration has the virtue 
of great simplicity, but does not allow for the ex- 
istence of a Poynting flux. They obtained fast ax- 
ial outflows but the luminosities required to drive 
them were as high as 10 5 — 10 6 L@. In a subse- 
quent work (Lery et al. 1999), this constraint has 
been relaxed. Only the necessary proportionality 
condition on the poloidal components of the ve- 
locity and magnetic field (v p cx B p ) was satisfied 
and the importance of the Poynting flux taken into 
account. The corresponding solutions were found 
to be faster and required a comparatively smaller 
source luminosity to be driven. In this work, we 
use the formalism and conditions from Lery et al. 
(1999) based on Poynting flux driving mechanism 
( y p{9) y<p{@)) t na t breaks the collinearity be- 
tween B and v. 

2.2.3. Self- similar radiation treatment 

In the early stages of star formation, most of the 
luminosity comes from the accretion shock created 
by the infalling material. Later on, the accretion 
rate reduces and the radiation is dominated by the 
protostar luminosity. Here, we use a simplified de- 
scription for radiation (Fiege & Henriksen 1996) 
that allows us to use an analytical self-similar ap- 
proach. 

Radiative diffusion: The steady-state energy 
conservation equation, that includes the mechani- 
cal and radiative energy fluxes, simply reads 



pw ( h + h rad + — + $ grav 



+ F 



rad 



= o 

(12) 
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where: 

• h is the specific enthalpy: h = u + p/p, with 
u the internal specific energy and p and p 
respectively the pressure and density. 

• h ra d is the specific radiative enthalpy. For an 
isotropic radiation field, the radiative pres- 
sure is p ra d = Urad/i, with u ra d the radia- 
tive energy density. Then, h ra d is defined as 

Kad = {Urad + Prad) I P = ^Prad/ P 

• &grav is the gravitational potential created 
by the central protostar. 

• Frad is the radiative flux. 

For the radiation independently, the radiative en- 
ergy conservation reads (Mihalas & Klein 1982) 



V • (Frad + Pvh 

r 



Kp 

= V • F raa 

C 



where c is the speed of light and k the opacity. It 
is re- written at the zevo-th order in (v/c) 



V-F r 



. 



(13) 



This equation decouples the mechanical energy 
fluxes from the radiative ones in Eq.(12). Fur- 
thermore, when taking the zcio-th order in (v/c) 
in Eq.(12), the h ra d term disappears. This means 
that the radiative pressure is neglected, i.e. the 
photons transfer no momentum to the gas. 

The previous statements are valid for any par- 
ticular form of the radiative flux F ra d- In this 
work, we use the radiative diffusion approxima- 
tion 



F r ad ^Prad 

up 



(14) 



where the radiative flux is directly linked to the 
energy density (then to the radiative pressure). 
This assumes that the mean free path of a pho- 
ton is very short compared to the characteristic 
length of the system. We use the black-body ra- 
diative pressure for p ra d which, at the temperature 
T is given by 

Prad^^-T 4 (15) 

3 c 

with a the Stefan-Boltzmann constant. 
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Fig. 1. — Opacity as a function of temperature 
computed with the code from Semenov et al. 
(2003) for a gas density of 10~ 18 g cm" 3 . 



Opacity: For the sake of simplicity, we choose 
an opacity defined by Kramer's law, 



K = Kq 



P 



1 g cm 



T 

Tk 



(16) 



where the values of the exponents a and b are de- 
fined according to the type of coolant considered. 
These parameters will be discussed in greater de- 
tail in section 4.3. For example, when the cool- 
ing is dust-dominated, a — and b = 2 (e.g. 
Pollack et al. (1985)). The constant k is deter- 
mined using a numerical code developed by Se- 
menov et al. (2003) which calculates opacities over 
a wide range of gas densities and temperatures, 
given a solar-type metallicity. This code includes 
different shapes and compositions for the dust par- 
ticles. In this work, we use the "homogeneous 
compact spherical dust" configuration to compute 
the opacity. Moreover, we focus on molecular out- 
flows which temperature range (10-100 K) is in the 
dust-dominated regime of opacity. In that case, 
the dust being such an efficient coolant, the opac- 
ity is independent of the density of the gas, i.e. 
a = 0. The opacities, obtained for a gas density 
of 10~ 18 g cm -3 , are plotted in Fig. 1 and can be 
fitted by a power-law. We found the Kramer's 
coefficient of the temperature to be b = 2 and 
k ~ 2 x 10~ 4 cm 2 g _1 for the dust dominated 
regime. 

Following Eq.(6), the self-similar radiative flux 
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reads 



Frad(r, 9) 



GM 
ro 



M 



ro 



f(6). (17) 



Using the self-similar expression of the tempera- 
ture T in Eq. (15) and combining Eq. (14) and 
(15), one finds a f to be a f = 6-3 + 2(a+l)(l/4- 
a). 

Furthermore, the inclusion of radiative diffusion 
in the model constrains the expression of the fidu- 
cial scale r . Indeed, the RHS and LHS multi- 
plicative constants of the self-similar expansion of 
Eq. (14) have to be equal: this leads to the follow- 
ing expression of r , 



To = 



3/c 
4a., c 



M 2+a (GM) 



4_§ / nm H \ 
V k B ) 



b — 4\ 5 + 6a + 2b 



(18) 

4<t/c with a the Stcfan-Boltzmann 



where a s 
constant. 

It is worthwhile to note that when diffusive ra- 
diation is discarded 1 , the fiducial scale is not ex- 
plicitly defined. In that case, the fiducial radius re- 
mains a free parameter and a characteristic length 
has to be chosen on observational or physical cri- 
teria (see Lery et al. (1999) and Sect. 2.3.2 for the 
discussion on this issue). 

2.3. Methods 

2.3.1. Integration 

Assuming self-similarity, the PDE system of 
HD/MHD equations is transformed into an ODE 
system where all variables depend only on the 
poloidal angle 6, 



dZj(fl) 
d9 



Mi im ■ 



(19) 



In the MHD (HD) case, eight (six) coupled equa- 
tions constitute the system. We consider the prob- 
lem as an initial value problem (IVP) and use the 
values of the variables (u r , ug, u<j,, jj,, y p , y^, G and 
f$) as input parameters at a given angle. Note 
that y p and are not present in the HD case. In 



1 i.e., there is no 0-componcnt in the radiative flux, but only 
a radial one, cf. the "virial isothermal" case in Fiege & 
Henriksen (1996). 



practice, we start the integration close to the ro- 
tational axis, typically at 6o = 10~ 2 rad. We pro- 
vide the initial values near the axis, i.e. a positive 
u r for the outflow, ug negative for the circulation 
pattern and u<f, chosen to be positive. We look for 
solutions covering # = #oto# = 7r/2 — t as the 
polar axis and the equatorial plane are excluded 
by the self-similar treatment (Fiege & Henriksen 
1996). 

2.3.2. Dimensions and the fiducial scale tq 

For each solution, dimensional quantities can be 
calculated using Eq.(7) to (11) and Eq.(17). From 
these equations, for any physical quantity, only 
the mass of the central object M and the fiducial 
scale r are required to compute the dimensional 
quantities at a given radius R. Without radiative 
diffusion, ro remains a free parameter (see Lery 
et al. (1999) for more details). 

In Fig. 2, the fiducial scale (from Eq.(18)) is 
plotted as a function of the Kramer's coefficients 
a and 6 for a one solar mass protostar and k = 



2 x 10" 4 



cm g 




Fig. 2. — Variation of r with the Kramer's opac- 
ity parameters. The values are obtained for a one 
solar mass protostar. 



The variation of r with the Kramer's param- 
eters is quite stiff and r quickly drops from ~ 
100 AU to as a increases and/or b decreases from 
the dust-dominated case. Note that ro is a lower 
limiting radius for the validity of self-similarity. 

For physical reasons, r should also be limited 
to the region where the hypotheses of the model 
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still apply. For example, close to the protostar, 
the radiation pressure could have non-negligible 
dynamical effects preventing the model to be ap- 
plicable there. 

2.3.3. Selection criteria of the solutions 

Firstly, we select solutions showing outflowing 
motions near the axis region and infall near the 
equatorial plane. This means that radial velocity 
u r must change sign once for an angle 6 opcn which 
is the opening angle of the outflow: for 9 > 9 open , 
u r is negative and trace the infall; for 9 < 9 open , 
u r is positive and the gas is flowing outward. 

Secondly, we use the self-similar expressions of 
the physical quantities to go from the dimension- 
less variables to the dimensioned ones: we define 
the mass of the central object (we typically choose 
1 M ), an observational radius (~ 5 x 10 3 AU) and 
a type of coolant (values of a and b in the opacity). 

Finally, we discard unphysical solutions, i.e. so- 
lutions that do not fit in the range of observations. 

3. Hydrodynamical solutions 

Pure hydrodynamical (HD) solutions can apply 
in situations where magnetic fields are either very 
weak (e.g., primordial universe) or not dynami- 
cally dominant for the movement of the gas (early 
stages of the collapse). Besides, studying a purely 
HD problem allows us to distinguish the effects of 
the radiation and density gradients from those of 
the magnetic field (see Sect. 4). 

We find HD solutions with both infall and out- 
flows and they all show common properties: i) the 
outflow is very narrow and ii) always quite slow 
(< 10 km s _1 for a one solar mass central object) 
and iii) infall is occurring in the rest of the domain 
in a quasi-spherical way. 

For illustration, a typical hydrodynamical solu- 
tion is plotted in Fig. 3: the left panel corresponds 
to the density contours and velocity field of the 
solution in the poloidal plane, and the right panel 
to different quantities plotted (at a fixed distance 
of 5000 AU) as a function of the poloidal angle 
9. The central object was assumed to be a one 
solar mass protostar. The Kramer's opacity pa- 
rameters are those of a dust-dominated opacity 
(namely, a = and b = 2). 

The infall-outflowing pattern is traced by the 



radial velocity: a negative radial velocity corre- 
sponds to an infall motion whereas a positive one 
is the signature of an outflow. Note also that vg 
keeps a constant negative sign, implying that all 
the material is being redirected outward in this 
particular solution. This will not be the case for 
solution with "net infall" (see Sect. 4.2). The 
velocity characterizes the movement of the gas 
towards the rotational axis. Here, the maximum 
reached by vq is ~ 0.45 km s _1 : this is a small 
value that implies a weak "collimation" of the flow. 
For this solution, the outflow has a velocity of ap- 
proximativcly ~ 1 km s _1 and a ~ 5° opening 
angle. The system is also almost non-rotating as 
«0 is close to zero over the entire domain. The val- 
ues of the number density (upper-left panel) are in 
the order of 10 6 cm~ 3 at 5000 AU and shows a in- 
crease towards the axis of rotation. However, the 
density at the axis is only twice as large as the one 
near the equator, emphasizing a quasi-spherical in- 
fall. At a given radius, the temperature (lower 
right panel) is also almost constant (~ 100 K at 
5000 AU) which characterizes an isothermal infall. 

Qualitatively, this solution presents some com- 
mon features with observations. First we have 
an outflow. We also find that the fastest mate- 
rial was also the most collimated (Bachiller 1996; 
Richer et al. 2000). This is precisely the case here 
where the positive part of the radial velocity (out- 
flow) that decreases when the angle from the axis 
increases. However, physical quantities (velocity, 
temperature and density) do not match with ob- 
servations. In particular, the density is an order of 
magnitude too high to fall in the observed density 
range (10 4 — 10 5 cm -3 ), and also the velocity is 
too small. 

In conclusion, we have shown that our models 
can produce outflows that can be launched from 
thermodynamical effects only. We have a sim- 
ple heated quadrupolar model for infall and out- 
flow. The qualitative characteristics of the hydro- 
dynamical solutions are similar to the observed 
protostellar outflows but are quantitatively too 
dense, too slow and too narrow for typical solar 
mass objects. 

4. MHD solutions 

Typical molecular outflows of low mass proto- 
stars have observed velocities ~ 20 km s _1 and 
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Fig. 3. — Left panel: a typical solution for the HD case represented in the poloidal plane {r,8). The iso- 
density contours and velocity field are plotted. Right panel: corresponding values of the velocity, density, 
temperature and radial radiative flux at a distance of 5000 AU from the central protostar, with respect to 
the poloidal angle 9. For this particular solution, the input values of the dimensionless quantities are: the 
velocities u r = 2.54, ug = —4.46 x 10~ 3 , = 9.05 x 10~ 3 , the density fi — 3.62 x 10~ 3 , the temperature 
= 1.81 and the 8— component of the radiative flux fg — 7.92 x 10 _1 . 



have a large range of initial opening angles: from 
< 30° for class to > 90° for class 1 objects 
(Bachiller 1996; Bachiller & Tafalla 1999). Their 
typical densities lie around 10 4-5 cm~ 3 . Our hy- 
drodynamical model cannot power such outflows 
and we now include the magnetic field in the 
model. 

There are two main issues that we wish to ad- 
dress with our magnetized model: 

• Class low mass protostars: they show pow- 
erful, highly collimated molecular outflows 
despite their very young age. 

• The formation of massive stars: they present 
very massive outflows (mass loss rates ~ 5 x 
10~ 3 Mq yr _1 ) that can be more massive 
than the central object itself and hardly be 
explained by the conventional jet- or wind- 
driven models (Churchwell 2000). 

4.1. Dust case a typical solution 

A typical MHD solution with pure transit is 
shown in Fig. 4. The left panel corresponds to the 



density contours and velocity field in the poloidal 
plane and the right panel to different quantities 
plotted at a fixed radius of 5000 AU as a function 
of the angle 6. 

As previously, the radial velocity traces the in- 
falling/outflowing motion of the gas. For this so- 
lution, the opening angle is 6> ope n = 1.1 rad « 60°, 
where v r changes sign. However, between 1.1 and 
0.3, the radial velocity is positive but very small 
and it is only beyond 0.3 rad w 17° that a non 
negligible outflowing motion is occurring. 

The outflow velocity is typically 10-20 km s _1 
for the most collimated part and the number den- 
sity varies from 10 5 to 4 x 10 5 cm~ 3 between 
< 6 < 7r/2. These values fit in the range of typ- 
ical values of molecular outflows (Bachiller 1996). 
The turning angle is the preferential angle for the 
toroidal magnetic field (middle right panel), re- 
sponsible of the collimation of the flow via the 
j z x B$ component of the Lorentz force. In conse- 
quence, vg, which traces the movement of the gas 
towards the axis, is also maximised at the turning 
point. Naturally, this correlation was not found 
in the HD case where the opening angle was very 
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Fig. 4. — Same as Fig. 3 but for a MHD solution with pure transit. For this particular solution, the 



initial values of the dimensionless variables are: u r = 1.45 x 10 , ug = - 
y p = 5.12 x 10 5 , y = 5.81 x KT 1 , fi = 5.93 x 10" 6 , 6 = 2.13 x 10" 1 and /< 
and the opacity parameters correspond to the dust case a = and 6 = 2. 



8.12 x 10~ 2 , m = 7.49 x 10~ 3 , 
= 1.353. The value of a is -0.22 



small but vg was maximum around 55 °. Remem- 
bering that i?^ oc in the model, v,p behaves 
identically to the toroidal field. 

The density (upper right panel) presents a 
strong poloidal gradient (dg) near the rotational 
axis, followed by an almost constant behavior for 
intermediate values of 9 and finally increases again 
when approaching the equatorial region. The 
strong gradient near the axis contributes to the 
acceleration of the outflowing material in this re- 
gion. This behaviour is not present in the HD 
solutions. The density structure of the present 
MHD solution appears oblate with respect to the 
rotational axis whereas it looks prolate for the 
hydrodynamical model. Such a geometrical dif- 
ference may be observationally distinguishable at 
high resolution. Indeed Motte & Andre (2001) 
have investigated the morphologies of Class and 
1 protostellar envelopes and have found that sev- 
eral sources (e.g. L1527, B335) have elliptical (or 
even more complex) density structure. They in- 
voke the presence of bipolar outflows to possibly 
be the reason of these asymmetries. This could be 
the case in our model. 

As in the hydrodynamical case, the tempera- 
ture is almost constant over the domain. However, 



the radiative flux plotted in the poloidal plane in 
Fig. 5, right panel, shows a strong anisotropy and 
is the highest in the axial zone and also contributes 
to the acceleration of the outflowing gas. On the 
other hand, the flux appears smaller in the equa- 
torial region than elsewhere, emphasizing that the 
equator is a privileged region for the accretion of 
matter onto the central protostar. That is not the 
case for the HD model (Fig. 5, left panel) where 
the radiation is almost spherical. 

The typical values for two models, hydrody- 
namical and magnetized, are gathered in Tab 1 
which summerizes the previous points. 

Knowing the velocity and the density for every 
angle 6 at a fixed radius R, it is possible to calcu- 
late mass rates. In particular, the infall mass rate 
is given by 

,tt/2 

M in = 2 x 2ttR 2 / v r (R, 9) p{R,9) sin 6AB 

J ^open 

(20) 

where v r (R,9) and p(R,9) takes the form of 
Eq.(10) and Eq.(7). The first factor 2 is present 
to take into the contribution to the infall from the 
other side of the equatorial plane. 

We compute the infall rates for different cen- 
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Table 1: Comparison between typical values of the HD and MHD solutions at a distance of 5000 AU from 
the central protostar. 




R (AU) R (AU) 

Fig. 5. — The radiative flux contours of typical transit solutions in the HD (left) and MHD (right) case 
between and 10 4 AU. At 5000 AU, the typical values are respectively 6 x 10 3 and 3 x 10 3 erg s _1 cm~ 2 . 



tral masses and plot the result in Fig. 6. Three 
solutions are treated: the one detailed previously 
in Fig. 4 (solid line) and two extreme cases (a low 
and a high density solutions, in dotted lines). The 
low and high density solutions are extreme in the 
sense that they are at the limit of the observational 
range but not a limit of the model (which can also 
produce denser or lighter solutions). Then, for a 
given central mass, our model predict infall rates 
located between the two dotted lines. 

First of all, the infall rate increases with the 
mass of the protostar and we get a typical value of 
M in w 3 x 10~ 6 M© yr -1 for a 0.1 M@ central ob- 
ject, which is in agreement with rates inferred from 
observations of low mass YSO (Bontemps et al. 
1996). It is important to mention here that the 
increase of the central mass does not correspond 
to the temporal evolution of the mass of the pro- 
tostar: we are studying a self-similar steady-state 
problem where the central mass is a dimensional 



parameter. Observationally, the higher the cen- 
tral mass is, the later in the pre-stellar evolution 
the protostar is, and as a consequence, the lower 
its accretion rate is. Here we are just studying 
the influence of the central mass on a given model 
and the increase of the infall rate is a consequence 
of the increase of the gravitational force from the 
source. 

The infall rate scales with M following a power- 
law. For this example, we have M- ln oc M 81 . This 
behavior comes from the self-similar form we use 
for the variables and, injecting Eq.(7) and Eq.(10) 
into Eq.(20) and using Eq.(18) for tq , one finds 

3 2(a + b-l/2) /- /o . n \ 

M in CX M$M~ 5+6a +a V (5/2+2a) _ (21) 

If a = and 6 = 2 (dust case), it simply reads 
M in oc iK 1 ""' and with a = -0.22 we find a 
slope of 0.81. In our model the smallest value a 
can tend to is —1/2 (a — —1/2 corresponds to 
pure radial infall) which would give a maximum 
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from solution of Fig. 4 

from 2 "extreme" solutions (high and low density) 



M (solar masses) 

Fig. 6. — Influence of the central mass on the in- 
fall rate for a MHD transit solution. The solid 
line corresponds to the specific solution presented 
previously. The dotted lines correspond to two 
"extreme" solutions of high and low density. 



slope of 1, Mi„ ocMin the dust dominated regime. 

With the purely transit solutions like the pre- 
vious one, the infall rate we calculate corresponds 
to the rate of gas moving towards the protostar. 
However, there is no net infall onto the central 
object in the sense that for these models, all the 
matter that falls is deflected in the outflow. This 
can be shown by writing the continuity equation 
Eq.(l) in the dimcnsionless variables and integrat- 
ing over all space: 



ir/2- 
0+e 



(1 + 2a)fj,u r H ^- (fine sin 6) sin9d9 

sin 9 dO 



(22) 



The first term corresponds to the total dimension- 
less mass rate over the entire domain. With our 
boundary conditions u$(0) = ug(7r/2) = 0, and fi 
finite, the integration of the second term equals 
zero. As a consequence, the total mass rate in the 
model is also equal to zero and all the infalling 
matter in diverted outward. In order to compare 
to observations, it is possible to compute the out- 
flow to infall rate ratio, / = M out /M in . For Class 
and 1 objects, the ratio is found to be smaller 
than unity typically, / ~ 0.1 — 0.3 (Bontemps et al. 
1996; Richer et al. 2000). In our case, for pure 



transit models, / = 1 by construction since all the 
material that is coming in has to be rejected in the 
outflow. 

Nevertheless, the inclusion of magnetic field 
greatly improves the properties of the solutions 
which now compare well with observational quan- 
tities, with fast, collimatcd outflows. 



4.2. 



Dust case 
infall 



a typical solution with net 



We now present solutions with a combination of 
transit and pure infall solutions similarly to Lery 
et al. (2002). We briefly review the interest and 
the properties of such solutions. 

The method of integration is the same as previ- 
ously except that, here, the domain is separated in 
two zones. The separation occurs at the angle 9 S . 
From #o to 9 S — e, we search for a transit solution, 
qualitatively identical to the one of the previous 
section. From 9 S + e to ir/2 — e, we look for a so- 
lution that presents the characteristic of an infall 
(v r < 0) but that is directed toward the equator 
and not toward the axis so that the gas is not de- 
viated. This condition is fulfilled if the 9— velocity 
is positive. At the end of the transit region, in 9 S , 
v$ tends to zero and is negative (requirement for 
transit). To integrate the system in the net infall 
region, we change the sign of v$, in order to get 
a strictly infalling pattern and use the end values 
of the transit solution, as initial conditions for the 
solution between 9 S and tt/2. 

There is a small region around 9 S which is not 
treated: 9 = 9 S is a location that cannot be 
reached (as for the equator in the pure transit 
solutions) and the integration numerically failed 
before reaching it. It is also important to men- 
tion that the two solutions, on both sides of 9 S 
are mathematically independent. The continuity 
set between the two zones is in consequence only 
apparent. 

A solution with the characteristics described 
above is shown in Fig. 7. As previously, we plot the 
density contours and the velocity field (left panel) 
as well as the variation of the physical quantities 
with 9 at a given radius (right panel) . We are con- 
sidering the standard case for which the opacity is 
dust-dominated. 

The two zones can be distinguished: i) the tran- 
sit zone, characterized by the change of the sign 
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Fig. 7. — Same as Fig. 4 but for a solution with net infall. 



of the radial velocity and by a negative theta- 
velocity. (i.e. from 9q — 10~ 2 to 9 S ~ 0.8 rad), ii) 
a net infalling region where v r < and > in 
the remaining of the domain. In this region, the 
gas is falling onto the equatorial plane so it can 
eventually accrete onto the central protostar. The 
untreated region between the two zones has been 
shadowed. 

Fig. 7 shows a solution with an opening angle 
of 9 opcn ~ 30° but solutions with smaller opening 
angles have also been found. The radial velocity 
for the outflow is 20-50 km s _1 and the average 
density around 2 x 10 5 cm~ 3 . Again, the solution 
appears almost spherical: each physical quantity 
does not show strong gradient except close to the 
boundary angles. Identically with the transit mag- 
netized solution, the minimum density and maxi- 
mum radiative flux is in the axis region: this allows 
the outward acceleration of the gas in this region 
and the radial velocity is then also maximum at 
the axis. 

The 9— velocity tends to zero at the axis and 
towards 9 S . These are the boundary conditions 
we require for transit solution and have in conse- 
quence to establish a pressure balance between the 
transit and net infalling regions. Near the equa- 
tor, vg does not tend to zero. This implies that the 
equatorial zone [tt/2 — e,7r/2] has to be seen as a 
sink that absorbs anything entering it. Looking at 



Eq.(22) for the net infalling zone (i.e. integrating 
between 9 S and 7r/2), one sees that a non-zero vq 
at one of the boundary (tt/2 in the present case) 
results in a non-zero global mass rate in the con- 
sidered region. This justifies a posteriori terming 
these solutions "net infall" models. For this cate- 
gory of models, the ratio / is then smaller than 
unity. For the particular solution shown here, 
M out - 2 x 10- 5 M yr- 1 and M ia ~ 2M out re- 
sulting in / = M out /Min ~ 0.5. This is a typical 
value for net infall models, but larger and smaller 
values can also be obtained by setting the net in- 
falling region as a smaller or larger fraction of the 
full model. 

4.3. Influence of opacity 

The solutions presented previously had a dust- 
dominated opacity. This is the situation in the 
present-day universe which is enriched in heavy 
elements, molecules and grains thanks to many 
generations of stars and stellar nucleosynthesis. 
However, in the primordial universe, star forma- 
tion occurred in an almost metal-free environment 
where the chemical composition was a mixture of 
hydrogen (H), helium (He) and very small traces 
of deuterium (D) and light elements (Li, Be, B). 
To understand the physics in the early universe, 
it is necessary to have accurate values of its opac- 
ities. Many studies undertaken on primordial star 
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formation used only the species involved in the for- 
mation of molecular hydrogen (and sometimes he- 
lium) (Stahler 1986; Omukai & Nishi 1998; Naka- 
mura & Umemura 1999) as H2 is the most abun- 
dant molecular species formed under primordial 
condition and is also the most effective coolant 
of the gas in the low temperature regime of star 
formation. However, HD molecules can also con- 
tribute to the cooling of the gas when H 2 be- 
comes inefficient (Palla 1999) and in a recent pa- 
per, Mayer & Duschl (2005) show that the inclu- 
sion of Li in primordial opacity calculation lead to 
significant changes, up to two orders of magnitude 
for T < 4000 K and conclude that the influence 
of Li on the different stages of population III star 
formation and evolution could be assessed. Here, 
we wish to study the behavior of our models when 
the optical properties of the fluid are those of a 
dust-free environment. 

The use of the Kramer's law form for the opac- 
ity prevents an explicit treatment of the cooling 
that includes all the species previously mentioned. 
However, the Kramer's opacity coefficients (a, b) 
can be interpreted in term of the physical process 
enabling the cooling. The emission rate per unit of 
volume can be written as e — AapnT 4 oc p a+1 T b+4 
(Boily & Lyndcn-Bell 1995). 

• a = corresponds to situations where a sin- 
gle efficient component (dust or CO) dom- 
inates the cooling, assuming its fraction by 
mass remains constant during evolution. 

• On the other hand, a = 1 gives the density 
dependence of the opacity when the cool- 
ing is due to two-particle processes and no 
coolant is dominant. For such a density de- 
pendence, the system is expected to be di- 
lute and cold so that the excitation levels are 
almost empty. 

• b measures the difference between a black 
body emission (oc T 4 ) and the one consid- 
ered here. In the case of dust, b = 2 as seen 
before. This appears to be an upper limit 
and b then decreases for other types of cool- 
ing, in particular b £ [—5/2, —1] for molecu- 
lar cooling in interstellar clouds (Goldsmith 
& Langer 1978). 

In order to study the influence of the opacity on 
the solutions of our model, we start from the dust 



case and then increase a from to 1 and decrease 
b from 2 to 1 so as to leave the dust-dominated 
regime. If no dust is present in the star forma- 
tion environment, then the main source of opac- 
ity would come from molecular line cooling . As 
mentioned previously, that case would see b in the 
range [—5/2,-1]. However, in this work, b does 
not take any value smaller than unity. This is due 
to the stiffness of the ODE system to solve. In 
fact, to study the influence of opacity on the so- 
lutions, we do not change any other parameters 
but a and b; for a given set of input parameters, 
a valid solution (i.e. filling the whole space and 
presenting circulation features) appears to remain 
valid only for b varying between 2 and 1. Below 1, 
the integration fails before 6 — tt/2 and the other 
input parameters should be slightly adjusted. Do- 
ing so would make it impossible to distinguish the 
effects due to opacity from those coming from the 
new input values. The limit b = 1 is not strict, we 
found solutions that could still be integrated with 
b < 1 and others that would fail earlier. 

By changing a and b progressively, we depart 
from the dust-dominated regime. It is however 
important to note that we cannot say what type 
of coolant it is we are modelling except that it 
is less efficient than dust and results in a smaller 
opacity. 

The study is made on pure transit solutions 
only and not on net infall ones as the disconti- 
nuity between the two regions of the latter was 
more difficult to handle. 

The results are shown in Fig. 8 and Fig. 9. 
The infall rate is calculated in the same way as 
in Sect. 4.1 and is plotted in Fig. 8 as a func- 
tion of the Kramer's opacity parameter a (b stays 
equal to 2). The infall rate increases as one leaves 
the dust opacity regime. Very simply, when the 
opacity decreases and the medium becomes more 
transparent, the photons escape more easily and 
are less efficient in slowing the infalling gas. 

In Fig. 9, the left panel represents the density 
for different sets (a, b) . The radial radiative fluxes 
are plotted on the right panel for the same sets of 
parameters. Furthermore, the quantities plotted 
are dimensionless (i.e., /z(0) and f r (9)) because as 
a and b are modified, the fiducial scale tq changes 
(see Eq.(18)): the influence of the opacity on the 
morphology of the solutions would be more diffi- 
cult to visualize as they would not start from the 
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Fig. 8. — Influence of Kramer's opacity parameter 
a on the infall rate. 
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Fig. 9. — Influence of the opacity parameters a 
and b on the dimensionless density and radial ra- 
diative flux. 



same point. 

The dimensionless density increases in the 
equatorial region as one leaves the dust-dominated 
regime. As it is the dimensionless density, this 
means that as a increases and b decreases, the ra- 
tio between the density at the equator and at the 
rotational axis increases. The radial radiative flux 
shows the opposite behavior: the flux is relatively 
smaller at the equator than in the axial region. 

It is not possible to generalize this behavior as 
not all the solutions that were found behave like 
that. For example, the morphology of the solu- 
tion presented in Sect. 4.1 was almost unaffected 
by the change of opacity. The only difference was 
observed for its dimensional quantities that would 
scale differently and, as a result, give a higher in- 
fall rate when dust is not the dominant coolant 
(see Fig. 8). However, it is of interest that such 
solutions exist: the higher density in the equato- 
rial region along with the weaker radiative field in 
this region could lead to the accumulation of mat- 
ter in the equatorial plane (i.e., to a more massive 
accretion disk) and may ultimately lead to the for- 
mation of more massive star for dust-free opacity 
regimes (such as the ones of the primordial uni- 
verse). 



5. General properties of MHD models 

The main feature of the model is the production 
of a heated pressure-driven outflow with magneto- 
centrifugal acceleration and collimation. An evac- 
uated region exists near the axis of rotation where 
the outflow is produced and the density system- 
atically increases with the angle from the axis. 
The streamlines passing close to the rotational 
axis are also the ones closest to the central ob- 
ject (see Fig. 11, upper left panel). The gas on 
these streamlines transits deeper in the gravita- 
tional well and is also the most vigorously heated. 
As a consequence, the radial velocity of the outflow 
is always maximum near the axis and decreases 
with increasing poloidal angle 9. 

5.1. Two families of solutions 

All the solutions that we have found have the 
same qualitative properties, except around their 
turning point. In this region, two distinct fami- 
lies of solutions appear as shown in Fig. 10. The 
rotational velocity (left panel) and the radial 
magnetic field B r (right panel) are plotted as a 
function of the poloidal angle 9 for different solu- 
tions. 

The solutions fall into one of the two cate- 
gories: i) slow rotators where the rotational ve- 
locity is small and almost constant over the entire 
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Fig. 10. — The rotational velocity (left) and radial 
magnetic field (right) plotted as a function of 9. 
The continuous and dotted lines correspond to the 
two types of solutions. For clarity we only show 
solutions with almost the same turning angle, but 
others were found. 



domain (dotted line) and ii) fast rotators where 
«0 increases and decreases with a strong gradi- 
ent around the turning point of the solution (solid 
line). At its maximum, peaks at about four 
times the v<f, value of the other solutions. Because 
we have oc v^,, the solutions that peak at the 
turning point also have a much stronger toroidal 
field there and as a consequence a larger pinch 
force towards the rotational axis. 

The radial magnetic field (Fig. 10, right panel) 
changes sign at the turning point of the solution 
(B r oc v r ) and it appears that the fast rotating 
solutions (solid lines) have a substantially higher 
radial magnetic field on the overall domain than 
the slow rotating models (from 2 to 16 times higher 
for the solutions shown). 

The other variables (not shown here: density, 
temperature, etc..) do not seem to be affected 
and no distinction can be made between the two 
types of solutions when looking at these quantities. 
No morphological differences appear between the 
two types of models and the density profiles all 
indicate the characteristic oblate shape of magne- 
tized solutions. Rotation and magnetic field seem 
coupled and the solutions with reasonable dimen- 
sioned quantities are in one of the two following 



categories: i) slow rotating and poorly magnetized 
or ii) faster rotating and highly magnetized. 

5.2. Energetics of the solutions 

In order to compare and contrast the two fami- 
lies of solutions, the evolution of the energy along 
a streamline is considered. We call Family 1/Fam- 
ily 2 the slow/fast rotating weakly/strongly mag- 
netized family of solutions found in the previous 
section. 

Streamlines follow dr/(rd9) = u r /ug in the 
poloidal plane. It is convenient to look at the lo- 
cal specific energy components normalized to the 
local specific gravitational energy (to have a di- 
mensionless quantity) and to define their sum as 

Be{r, 9) = —J— ( E k + E grav + V - + E magn ) 

E grav V P J 

(23) 

where Ek is the specific kinetic energy, E glav = 
—GM/r the specific gravitational potential en- 
ergy, and -E m agn = B 2 /8irp the specific magnetic 
energy. Using the self-similar form for the quanti- 
ties, Be becomes 

Be W = I K + „ S + »y-l+e+^(f + |) 

(24) 

Note that in the hydrodynamical case, the Bernoulli 
theorem implies that Be is constant along a 
streamline. 

In Fig. 11, upper left panel, streamlines are 
plotted for two solutions, one from each family. 
The central object is located at the origin. It is 
interesting to note that for streamlines integrated 
from the same point (located in the outflowing re- 
gion), the solution from Family 2 transits closer 
to the central object, deeper in the gravitational 
well than the one from Family 1. This is a general 
trend that has been observed for many solutions. 
The three other panels in Fig. 11 show energy com- 
ponents plotted along these two streamlines. 

On lower left and right panels the different spe- 
cific energy components as a function of the po- 
sition along the streamline of the solution from 
Family 1 and 2 are plotted. The origin of the 
curvilinear abscissa s is chosen at the "infalling 
end" of the streamline. 

Far from the turning point, the main contri- 
bution to the kinetic energy (E^) comes from the 
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Fig. 11. — Family 1 and 2 refer respectively to the 
slow and fast rotating solutions found previously. 
Upper left: poloidal streamlines. Upper right: ra- 
tio of the kinetic to the magnetic energy as a func- 
tion of the curvilinear abscissa s. s — is chosen 
at the infalling end of the streamline. Lower left: 
energetics of a solution along a streamline for a 
solution from Family 1. The specific energy com- 
ponents have been normalized to the gravitational 
potential energy. Lower right: same as lower left 
for a solution from Family 2. 



radial velocity v r (vq and have much smaller 
values). The absolute value of v r rapidly decreases 
when approaching the turning point (See Fig. 4) 
and so does Ek (dashed line), for both Family 1 
and 2 solutions. Near the turning point, how- 
ever, the contribution of ve,<j> to Ek becomes dom- 
inant: i) for the slow rotating solutions (Family 
1 - lower left panel), does not change much 
on the domain and Ek reaches a minimum at the 
turning point, ii) but for the fast rotators (lower 
right), ^0 quickly increases at the turning point 
and Ek reaches a local maximum. The magnetic 
energy qualitatively follows the same behavior as 
B r fi cx tv,i3 and cx v^. 

Finally, the ratio of the kinetic to the mag- 
netic energy is plotted on the upper right panel 
in Fig. 11, again for the two types of solutions. 
Qualitatively, they behave alike, and -Efc/-E mag n is 
in both case minimum at the turning point. This 
conforms with intuition as one expects a larger in- 
fluence of the magnetic field near the central pro- 



tostar. However, at the turning point, -E mag n is 5 
times smaller than Ek for the slow rotating solu- 
tion while it is comparable to Ek for the fast ro- 
tating model. For clarity only two solutions have 
been plotted here, but it has been found a general 
behaviour that (Ek / E magD )^™ lly2 ~ 1. Further- 
more, the same ratio for Family 1 solutions has 
been found many times much greater (5-100) than 
unity. 

In conclusion, the behaviour of the two types 
of solutions is similar far from their turning point, 
but then differs significantly: for the slow rotat- 
ing models, the magnetic energy is always much 
smaller than the kinetic energy whereas they be- 
come comparable for the fast-rotating solutions. 
The magnetic field appears then to be responsible 
for a closer transit around the center of the latter 
ones. Despite this difference between the two fam- 
ilies of solutions, the properties (velocity, density) 
of the outflows they power are similar. 

6. Discussion 

6.1. Behavior at large distance 

After the main early phase where infall domi- 
nates the dynamics, the protostar deposits linear 
and angular momentum and mechanical energy 
into its surroundings through its jets and molec- 
ular outflow. It is therefore important to discuss 
the behavior of the outflows at large distance. In 
the case of low mass stars, the accretion-ejection 
engine should dominate the global dynamics and 
produce relatively weak molecular outflows. In 
that case, the present model may apply only for 
a short stage at the earliest phases of low mass 
star formation. On the other hand, if the physics 
of the outflow formation remains the same regard- 
less the mass of the protostar, the dynamics of 
massive objects should be dominated by the tran- 
sit and should have strong and heavy molecular 
outflows, according to the present model. In that 
case, the transit model plays a major role in the 
formation of massive stars. Moreover, the tran- 
sit model brings an interesting new feature for the 
interaction between molecular outflows and jets. 
Indeed, a lot of gas is brought to the axial region 
that is not at rest and is already stratified when 
the atomic jet starts propagating away from the 
protostar. It will drastically change the jet propa- 
gation and stability properties. This demonstrate 
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the importance of understanding of the history of 
outflows in order to model correctly the behavior 
of jets. 

Also, according to our model, at larger distance, 
p oc r 2 «- 1 /2 ; B K r «-3/4 ) and y K r -i/2_ The 

magnetic field and the density decrease faster than 
the velocity. For example for a = —1/4, p oc r _1 , 
B oc r _1 , and u oc r~ 1//2 , and the outflow be- 
comes purely ballistic (Lery et al. 2002). Only 
large bow-shocks due to the long time-scale vari- 
ations of the source should survive. Ultimately, 
the combined effects of the reduction of the ion- 
ization and temperature in the giant bows, and of 
the lack of ambient medium, may make the out- 
flows fade out and finally become very difficult to 
observe on distances of 3 to 10 pc. However, the 
momentum that they may transfer to the ambi- 
ent medium could help to maintain turbulence on 
large scales. 

6.2. The outflow evolution 

Of our steady-state model, any temporal evolu- 
tion is beyond the scope. However, we may spec- 
ulate on the evolution of YSOs by considering the 
slow quasi-steady evolution of our model. With 
the new insight given by the model, we suggest 
some physical justifications for the early evolu- 
tionary sequence of both low and high mass YSOs 
which is usually defined on the basis of an empir- 
ical scheme only. 

During the Class stage, most of the mass of 
the system is still located in the infalling envelope. 
The source is totally embedded and the transit of 
matter around the center object can start during 
the formation of the stellar core. This gives rise to 
a fast, powerful and collimated molecular outflow. 
However during this early stage, the central object 
may not have already developed the jet that is 
commonly thought to drive the molecular outflow. 

When reaching the early Class 1 stage, the cen- 
tral object is surrounded by both a disk and a 
diffuse circumstcllar envelope. A jet eventually 
emerges from the newly formed inner accretion 
disk as the transit continues. The jet may then 
pressurize and entrain the material in the axial re- 
gion, pushing the transit flow away from the axis, 
thus reducing the molecular outflow velocity and 
increasing its diameter. Therefore, molecular out- 
flows will appear to be less powerful and less col- 



limated with time. 

For low mass protostars, in the latest stages of 
the pre-stellar evolution, the central object should 
not be embedded in an envelope and our model 
cannot apply. The molecular outflow continues to 
spread out and slow down and is eventually over- 
taken by the central jet still fed by the accretion 
disk. In the case of high mass stars, the tran- 
sit features may remain in action until the central 
star forms and reach the main sequence. Either 
the increasing opening angle of the molecular out- 
flow in the low mass case, or the ignition of the 
central object may ultimately stop the accretion 
by suppressing any contact between the reservoir 
formed by the molecular cloud and the accretion 
disk or the star. 

7. Summary 

In this paper, we have presented a MHD 
model that applies to infall and molecular out- 
flows around young stellar objects. The model 
is based on the radial self-similarity assumption 
applied to the basic equations of ideal, axisym- 
mctric and steady-state MHD, including Poynting 
flux. Instead of the usual mechanisms invoked for 
the origin of molecular outflows (underlying jet or 
wind), the outflow is powered by the infalling mat- 
ter through a heated quadrupolar transit pattern 
around the central object. 

Solutions without magnetic field have been 
found and they appear almost spherical (quan- 
tities do not vary much with the poloidal angle). 
Although they do not compare well with obser- 
vation ranges (too dense and too slow outflow), 
they indicate that thermodynamics is a sufficient 
engine to generate an outflow. They could apply 
to star formation in the primordial universe or 
in the very early stage of star formation where 
magnetic fields may not be dynamically impor- 
tant (however, in the latter case, the point mass 
gravity field should be replaced by self-gravity). 

The magnetized solutions show dynamically 
significant density gradients in the axial region, 
precisely where the radial velocity and collimation 
are the largest. Their radiative field is also highly 
anisotropic and the highest at the rotational axis. 
It contributes to the ejection of the gas. Quanti- 
tatively, the pure transit solutions compare well 
with observations of outflows: velocities of a few 
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tens of km s 1 and typical density of 10 4 5 cm 3 . 

However, for these models all the infalling gas 
is deviated into the outflow - there is no net in- 
fall. ft is possible to obtain, from the same set 
of equations, a purely infalling region around the 
equatorial plane. The gas in this region is not de- 
viated and will likely end up on the accretion disk 
around the protostar. When changing the scaling 
and putting a more massive central object, the in- 
fall rate also increases. This indicates that our 
model favors the large accretion rate scenario to 
form massive stars. 

The influence of the opacity on the transit so- 
lutions has been studied. When leaving the dust 
dominated regime, the fiducial scale of the model 
decreases which results in an increase of the infall 
rate. Furthermore, the morphology of some solu- 
tions is also affected: the density ratio between 
the equator and the axis increases while the radial 
radiative flux ratio decreases. As a consequence, 
when dust does not dominate the cooling, such as 
in the primordial universe, matter could be more 
easily accumulated in the equatorial region and 
form more massive stars . 

We suggest that molecular outflows are domi- 
nated by the global transit of material around the 
protostar (except for a thin layer surrounding the 
central jet - if present, where the dynamics is gov- 
erned by entrainment) and that the same process 
occurs from low to high mass forming stars. The 
present work also suggests that radiative heating 
and magnetic field may ultimately be the main 
energy sources driving outflows during star for- 
mation, at the expense of gravity and rotation. 

Finally, although the details of the outflow 
mechanisms may be peculiar to individual ob- 
jects we believe the infall-outflow circulation arises 
naturally given accretion, and thus could also be 
present in other astronomical objects such as ac- 
tive galaxies and around suitably placed compact 
objects, such as neutron stars. 
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